Nature Computational Science

1Clinical metagenomic sequencing is increasingly used to detect pathogens directly from clinical samples, but sparse and uneven genome recovery caused by low microbial biomass and stochastic fragment loss makes species-level validation difficult. Under these conditions, the same observed coverage may arise either from genuine genome-wide representation or from reads confined to conserved loci shared among related taxa. Current heuristic thresholds based on read counts, breadth or related alignment summaries target the realised coverage rather than the generating process. We show that, under hidden fragment loss, coverage is the wrong inferential object for uniformly valid species-level inference. What remains informative is the spatial dispersion of surviving reads across genomic coordinates. This change of inferential object allows us to ask how much of the candidate genome must be represented in any source capable of reproducing the observed evidence. This yields the Metagenomic Alignment Validation Index (MAVI), a conservative lower bound on the fraction of the candidate genome required to explain the observed read dispersion. MAVI therefore provides an error-controlled, interpretable measure of genomic confidence: the minimum genomic content a competing source would need to share with the candidate to mimic the observed evidence. Evaluated across 643 long-read metagenomic samples, including 280 clinical samples, MAVI reproduced independently supported organism calls, resolved sparse low-coverage ambiguities and exposed workflow-level collapse of stochastic dispersion, without relying on heuristic coverage thresholds.

Proteins are dynamic molecules existing in diverse conformational states underlying their biological functions. Although recent approaches have enabled diverse conformational sampling by emulating molecular dynamics simulations, perturbing evolutionary information, or steering internal mechanisms of structure prediction models, predicting conformations resulting from major domain motions or motions that occur over long timescales still remains a challenge. To this end, we introduce GNMCADS, a conformational sampling strategy that enhances the diversity of protein diffusion models by selectively annealing the conditioning signal guided by the intrinsic dynamical organization of the sampled protein. Further, we implement GNMCADS in the diffusion module of AlphaFold3, enabling the generation of diverse protein conformations. When benchmarked across 92 proteins that include 54 class A GPCRs, 15 transporters, and 23 proteins with major domain movements, GNMCADS exhibits improved sampling diversity compared to other current conformational sampling methods.

Biobank-scale genomic analyses are increasingly constrained by computational costs, as hundreds of thousands to millions of samples and variants must be analyzed together. The genotype representation graph (GRG) compactly encodes population genetic variation to accelerate computation, but the current approach does not exploit modern accelerator architectures. This work introduces Mikado, a new methodology for expressing GRG-based computation using sparse linear algebra primitives. Under a reverse topological ordering of the graph nodes, the GRG adjacency matrix is strictly block-lower-triangular, and the genotype matrix-vector product becomes a sparse triangular solve that can be further decomposed into a pipelined sequence of blocked sparse matrix-vector multiplies. By decoupling computation from graph representation, our approach exposes fine-grained parallelism and enables hardware-optimized sparse primitives on GPUs. Mikado achieves an order-of-magnitude speedup and cost savings for PCA and BOLT-LMM compared with the original GRG traversal approach, including on All of Us cohorts. It provides a scalable, hardware-portable, researcher-friendly tool for population genetics at biobank scale.

Population-scale sequencing now produces variant call sets with thousands of samples and millions of sites, making post-calling analysis a recurring bottleneck. Because the Variant Call Format stores every field of every record together, a tool answering a field-limited question still parses the unused annotations, FORMAT blocks, and per-sample values, which costs time without changing the result. We present VariantFlow, a command-line engine built on selective execution, decoding only the fields the current filter, statistic, or projection requires and preserving original records where possible. On correctness-matched benchmarks on the 1000 Genomes 3,202-sample high-coverage dataset, VariantFlow computed whole-chromosome missingness 8.0-8.5x faster than VCFtools at a constant 8.6 MB of memory, the chromosome 1 summary taking 128 s against 1031 s. An index-assisted FILTER=PASS predicate ran 123-273x faster than bcftools where block metadata excludes most of the file. Further gains, reported in Table 1, cover FORMAT-rich filtering on record subsets, linkage disequilibrium, and export-once Parquet queries served through DuckDB. Core population-genetic outputs were byte-identical to VCFtools, per-individual missingness and the site-frequency spectrum matched scikit-allel exactly, and the missing-data-aware pi and dxy estimators reproduced pixy's pairwise counts in every window. By decoding only the fields each operation requires, VariantFlow brings large-cohort post-calling analysis within reach of exploratory work on a single workstation. It complements rather than replaces bcftools, HTSlib, GATK, and VCFtools, and should be of value wherever large cohort VCFs are repeatedly summarised, filtered, or exported. VariantFlow is open source (Rust, MIT OR Apache-2.0) at https://github.com/ehsanestaji/VariantFlow; version 1.5.0 is archived at doi:10.5281/zenodo.21198172.

Medical AI models have made a great impact on biomedical research and real-world clinical applications, but conducting interdisciplinary medical AI research remains challenging, requiring close collaboration between clinicians and AI experts. Recent advances in large language models (LLMs) and autonomous code agents present an opportunity for low cost medical AI development, where clinicians can build AI tailored to their own research questions, even without continuous support from dedicated AI experts. However, enabling code agents to autonomously tackle complex multimodal medical AI development tasks requires clinicians to construct and supervise an AI research loop with detailed technical specifics, demanding substantial expertise in AI and computer science that they often lack. To address this challenge, we introduce the Medical AI Research Loop Agent (MARLA), an agentic framework that completely abstracts the construction and supervision of medical AI research loops from clinicians. Given a clinician-defined research intent, MARLA automatically translates high-level research goals into executable hierarchical research loops, decomposes them into verifiable sub-loops, and specifies the models, datasets, tools, and evaluation protocols required for each task. During execution, MARLA coordinates specialized code agents, monitors progress, diagnoses failures, and iteratively refines research strategies based on experimental feedback to drive the research process toward optimal outcomes. We evaluate MARLA on multimodal medical AI tasks that require closed-loop conversion from high-level clinical study objectives to trained and validated AI models. The results demonstrate MARLA's ability to autonomously conduct complex medical AI research while substantially reducing the need for AI expertise.

Accurate long-term forecasting of cognitive trajectories across the Alzheimer's disease continuum is essential for early intervention, personalized prognosis, patient stratification, and clinical trial enrichment. Despite the promising predictive performance of recent longitudinal forecasting methods, they remain largely data-driven, struggle with irregularly sampled, incomplete longitudinal data and often neglect established disease biology, leading to biologically implausible trajectories. To address this, we propose a biologically constrained continuous-time framework for long-horizon cognition forecasting from limited baseline observations. The proposed method models the complete amyloid-tau-vascular-neurodegeneration-cognition (ATVNC) cascade using hierarchical Neural ODEs with biologically motivated monotonicity constraints. Each pathological stream is governed by a dedicated Neural ODE initialized from irregular longitudinal observations using a GRU-D encoder, capturing intrinsic disease evolution while being modulated by directed upstream pathological influences. A bounded cognition readout ensures physiologically valid cognitive score (MoCA) predictions, while teacher-student knowledge distillation improves learning from sparse longitudinal supervision. Evaluated on the ADNI dataset, the proposed framework achieves a long-horizon extrapolation MAE of 2.06 on 188 held-out participants while eliminating biologically implausible trajectory violations. It further demonstrates robust zero-shot cross-cohort generalization on OASIS-3 (MAE 2.68 on 300 participants), with fine-tuning improving MAE to 1.90. The model also supports prognostic enrichment for Alzheimer's clinical trials, achieving up to 2.70x enrichment over the cohort base rate. These results demonstrate that embedding biological disease mechanisms within continuous-time deep learning improves the accuracy, biological plausibility, and clinical utility of long-horizon cognitive forecasting. The code is publicly available at: https://github.com/PonDeepika/BEACON.

Mitochondrial DNA (mtDNA) is increasingly recognized as an important contributor to human disease and population variation, yet most genomic biobanks do not provide standardized mtDNA variant datasets despite abundant mitochondrial sequencing reads in existing whole exome and whole-genome sequencing data. We developed a scalable framework based on the Mitoverse mtDNA Server 2 Fusion workflow to generate harmonized, analysis-ready mtDNA resources across diverse biobank infrastructures. The framework was implemented in the Mount Sinai Million Health Discoveries Program (54,151 participants) using the native Nextflow workflow and adapted for the All of Us Research Program (197,361 participants) using a custom cloud implementation that preserved the same analytical strategy. Across 251,512 participants, the framework generated standardized mtDNA datasets containing 12.9 million variant observations suitable for downstream genomic and electronic health record linked analyses. This framework enables reproducible, population-scale mitochondrial genomics across institutional and national biobanks without requiring additional sequencing or development of new variant calling methods.

Inferring the history of effective population size (Ne) from genomic data is fundamental to population genetics, but model-flexible methods have not kept pace with biobank-scale samples. We present Stairway Plot v3, which infers a nonparametric Ne (t) trajectory with bootstrap confidence intervals directly from the site frequency spectrum (SFS). By recasting the optimization as a convex problem and combining inferences across multiple sub-sampled projections, the method scales to 100,000 haploid sequences on commodity CPU hardware without a GPU. Larger samples sharpen the resolution of recent demography, and we apply the method to five human populations including 54,302 Japanese individuals, to our knowledge the largest sample used for nonparametric SFS-based Ne (t) inference.

Single-cell transcriptomic atlases resolve thousands of putative brain cell types, yet common analytical workflows often rely on low-dimensional embeddings that can distort neighborhood relationships and obscure continuous variation between closely related populations. We introduce single-cell Discrete Morse Graph Construction (scDMGC), a deterministic framework based on topological data analysis and discrete Morse theory that extracts a compact graph from neighborhood structure in high-dimensional gene-expression space. Because graph edges are derived from the original k-nearest-neighbor complex, scDMGC provides an interpretable multiscale representation of transcriptomic organization without requiring a low-dimensional embedding. Local maxima represent coherent transcriptional states, saddles quantify their separation, and gradient paths capture intermediate states and continuous transitions. Applied to cortical and hippocampal datasets, scDMGC recovers established inhibitory-neuron classes and elucidates transcriptional gradients. In mouse whole-cortex data, Morse graph structure evaluates cell-type hierarchy and type robustness. Persistence across scales provides a quantitative measure of cell-type identity and discrete versus continuous relationships. Finally, in Alzheimers disease single-nucleus data, scDMGC distinguishes changes in cell-type abundance from disease-associated shifts in transcriptional state. scDMGC therefore provides a multiscale framework for quantifying whether transcriptional populations form persistent cell types, how strongly external annotations are supported by intrinsic data structure, how cell types relate to one another, and where discrete organization transitions into continuous variation.

Short tandem repeats (STRs) are a rich and highly polymorphic source of human genetic variation, but representing and genotyping them in pangenome graphs remains challenging. Explicitly encoding each STR allele as a separate graph path results in increasingly complex local structures as cohort diversity increases, leading to larger index sizes and requiring significant resources for graph reconstruction when new alleles are introduced. Here, we propose STR-PG, a topologically decoupled genome-wide framework that separates stable locus representation from scalable STR allele content. STR-PG uses topologically fixed pointer nodes to represent each target locus, while allele sequences, repeat counts, motif annotations, and population frequency metadata are stored in an external registry. Short reads are mapped to STR loci via syncmer-based flanking anchors, and genotyping is performed within a locus-specific candidate space using allele-level alignment likelihood and Bayesian inference. Newly supported alleles can be integrated through registry-level updates without the need to rebuild the graph structure. Evaluations using simulated whole-genome sequencing data, 1000 Genomes Project (1kGP) samples, and r real whole-exome sequencing data from matched whole-blood-cell controls demonstrate that STR-PG maintains accurate genotyping results across various STR classes, reproduces expected population structures, and substantially reduces the computational cost of integrating additional alleles. STR-PG provides a compact and scalable framework for population-scale STR analysis using short-read sequencing.

Mendelian randomization (MR) uses genetic variants as instrumental variables for causal inference, but its validity is often undermined by horizontal pleiotropy and linkage disequilibrium (LD) among instruments, which current methods typically address by discarding correlated or outlying variants. Here we present ECC-MR, an MR framework that brings ideas from error-correcting codes into causal inference: a parity-check matrix constructed from LD structure places the causal signal in its null space, so that pleiotropic effects become a decodable sparse error syndrome rather than a nuisance. We establish theoretical guarantees for consistency and stable recovery under sparse pleiotropy, and show in simulations that ECC-MR outperforms existing methods in bias, efficiency and Type I error control, while flagging regimes where its assumptions fail. Applied to three exposure-outcome pairs spanning weak to strong pleiotropic burden, ECC-MR yields concordant causal estimates while retaining all correlated instruments, and its decoded pleiotropic variants recover known hotspots and drug targets, corroborated by colocalisation and phenome-wide association analyses. ECC-MR thus turns LD and pleiotropy from liabilities into resources for robust causal inference and mechanistic discovery.

Focal cortical dysplasia (FCD) is among the leading causes of drug-resistant focal epilepsy and can be challenging to detect on routine MRI due to the subtle and heterogeneous appearance of lesions. Motivated by the scarcity of clinically representative cohorts, we investigate the feasibility of simulating FCD-like lesions in healthy MRI using a composition of hand-crafted image transformations. We develop SynthFCD, a phenomenological FCD simulator, and create 2,164 synthetic examples from healthy controls in the heterogeneous FOMO300k dataset. We then pretrain deep learning (DL) segmentation models on synthetic lesions and fine-tune them on real FCD cases from the University Hospital Bonn (UHB) cohort. Comprehensive evaluations show consistent improvements across simulator presets, with detection-rate gains of up to 18.5 percentage points over training from scratch and 7.4 points over self-supervised pretraining, while achieving state-of-the-art segmentation performance (0.366 DSC) on the challenging held-out UHB benchmark. These findings suggest that simulation-based synthetic data generation can bridge existing radiological knowledge and advances in DL to provide meaningful supervisory signals that could ultimately aid the detection of underrepresented and diagnostically challenging cases.

Predicting the course of Alzheimer's disease for individual patients remains a major challenge due to the heterogeneity of disease expression and the sparsity of longitudinal data. We introduce a variational Disease Progression Score (DPS) framework that maps multimodal biomarker dynamics (Cerebrospinal fluid, neuroimaging, and cognitive assessments) onto a continuous latent timeline with quantified uncertainty. The framework combines a neural encoder, which infers subject-specific progression parameters from demographic and clinical features, with a cascade of logistic functions structured according to the amyloid cascade hypothesis. Applied to the Alzheimer's Disease Neuroimaging Initiative (ADNI) cohort, the inferred timeline separated diagnostic groups it never observed (AUC 0.98 for cognitively normal vs Alzheimer's Disease), and the estimated cascade strengths and biomarker orderings were consistent with the established sequence of Alzheimer's pathology. The model produces individualised prognoses for previously unseen subjects from baseline data alone, with 95\% credible intervals achieving 89-98\% empirical coverage across biomarkers, and these predictions can be dynamically refined as new observations become available. The framework thus provides a biologically interpretable, uncertainty-aware index of disease severity, offering a probabilistic foundation for patient-level prognosis and precision monitoring in Alzheimer's disease.

Brain rhythms organize neural activity in time. When a shared oscillation is broadcast across a population of neurons, the phase at which the neurons fire with respect to the oscillation can itself carry information, a strategy the brain appears to use for navigation, memory, and sensory perception. In such a phase-of-firing code, information is conveyed not by how rapidly a neuron fires but by when it fires relative to the phase of the ongoing rhythm. Models that implement such phase-of-firing codes consider the neurons as independent. How such codes arise in cortical and hippocampal circuits, where neurons are densely interconnected rather than independent, remains unclear. Here we show that recurrent connectivity does not merely permit phase-of-firing codes but actively determines which patterns the network can hold onto. Using networks of model neurons wired in a manner determined by the constraints of a Sudoku puzzle, we find that the network supports an enormous repertoire of distinct rhythmic firing patterns, each corresponding to a valid Sudoku solution and selectable by brief targeted inputs against a common background oscillation. Some patterns persist after the input is removed while others dissolve. We show that a factor deciding this selectivity is network symmetry. Patterns that respect the symmetries of the wiring are stable, while those that break them are not, and small perturbations can switch the network between stable patterns. Adjusting synaptic weights to remove asymmetry stabilizes patterns that were previously unstable, giving the network a tunable memory. More broadly, our results suggest that a common oscillatory drive, targeted biases to the neurons, and recurrent network structure, act together as complementary control variables. The oscillation sets up a temporal scaffold, the biases select spatiotemporal patterns, and the connectivity determines the stability of the chosen patterns. Tuning these complementary variables can predictably and flexibly reshape the repertoire of activity a circuit expresses. Significance StatementMany brain regions generate rhythmic activity, and the precise moment at which individual neurons fire within each rhythm carries information that the brain can harness to encode memory, perception and location as an animal navigates its environment. In such a phase-of-firing code, information is conveyed not by how rapidly a neuron fires but by when it fires relative to the phase of the ongoing network rhythm. Existing theories of how phase-of-firing codes arise have modeled neurons as independent, leaving open how such codes emerge in the densely interconnected circuits of cortex and hippocampus. Using model networks wired to reflect the constraints of the Sudoku puzzle, we show that recurrent connectivity supports a vast repertoire of distinct rhythmic firing patterns, and that the symmetries of the wiring determines which of those patterns the network can stably hold. This identifies network symmetry as a general principle relating circuit architecture to temporal coding.

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

MotivationSimulation frameworks for statistical and population genetics often separate reference-panel genotype generation, admixture modeling, phenotype simulation, and scalable genomic storage. This limits construction and evaluation of large synthetic cohorts with controlled donor composition, local ancestry, and phenotype architecture. ResultsWe developed FORGEPHAST, a Python framework for reference-panel-based synthetic haplotype and phenotype generation. FORGEPHAST supports homogeneous mosaic and pulse-admixture generation from user-defined donor panels, integrated genotype quality control, and phenotype simulation with population- and local-ancestry-specific genetic effects. Its HapStore backend provides haplotype-native, Zarr-based storage for scalable numerical analysis. We demonstrate preservation of key donor-panel population-genetic properties and use FORGEPHAST to compare polygenic risk score methods across 54 simulated phenotypic architectures at biobank scales, showing that relative performance depends on phenotype architecture, ancestry, and implementation-specific variant coverage. Availability and implementationFORGEPHAST is implemented in Python and is available as an open-source package at https://github.com/sriraj-lab/forgephast.

Generating diverse, all-atom protein conformations remains challenging, especially for recovering alternative functional states and modeling conformational ensembles. We present trFlow, an ultrafast approach for generating all-atom conformations using geometry-guided flow matching. trFlow leverages three core innovations: (1) a hybrid single- and multi-step sampling scheme that balances accuracy and speed; (2) a decoupled representation-generation architecture that maximizes computational efficiency; and (3) fine-tuning on NMR ensembles to better capture conformational heterogeneity. Extensive benchmarking across dual-conformation proteins, NMR ensembles, and a CASP16 dynamic target shows that trFlow outperforms competing methods in recovering challenging states and reproducing experimental heterogeneity, while being an order of magnitude faster than existing all-atom approaches. These results establish trFlow as a practical and scalable framework for large-scale generation of all-atom conformational ensembles.

Linkage disequilibrium (LD) information from one's own dataset is considered optimal for downstream analyses such as statistical fine-mapping. However, the computational complexity ({approx}N2 / 2 computations for N variants) lead most studies to use external reference panels, such as from 1000 Genomes. To capture LD across all variant pairs in biobank-scale whole-genome sequencing (WGS) datasets with hundreds of millions of variants, new computational strategies are essential. We present a novel approach that uses multi-GPU distributed computing to compute R2 for every variant pair in a dataset. On chromosome 22 (1.8 million variants) of the 30x WGS 1000 Genomes dataset, our method using eight consumer-level 12GB GPUs (NVIDIA RTX 2080TIs) takes <10 minutes, while the same calculation with PLINK using a high-end 64-threaded CPU (Intel Xeon Gold 6338) takes >80 minutes, corresponding to a {approx}8x speedup. We further demonstrate true biobank-scale performance in the Penn Medicine Biobank (PMBB; 57,170 samples), computing chromosome 1 LD ({approx}1.38 million variants) in 49 minutes versus 22.7 hours for PLINK, a {approx}28x speedup. With this method, we successfully computed on the aforementioned 30x WGS 1000 Genomes dataset ({approx}120 million variants and {approx}2500 samples) the entire LD matrix (>1e16, or 10 quadrillion elements) in under 6 hours using 512 NVIDIA 40GB A100 GPUs on the Department of Energy Argonne Leadership Computing Facility Polaris Supercomputer. We make this tool, coded in Python using CuPy, publicly available. Using this tool, researchers can leverage the full extent of their genomic data without relying on external LD reference panels and acquire more accurate, population-specific findings, particularly for groups underrepresented in existing databases.

Accurate prediction of antigenic distance between influenza A/H3N2 strains is essential for timely vaccine strain selection, yet traditional hemagglutination inhibition (HI) assays are labour-intensive and limited in throughput. We present FluQRC, a hybrid Quantum-Classical Reservoir Computing framework for sequence-based antigenic distance prediction. FluQRC integrates three novel components: (1) a differentiable gated property ranking network for data-driven property selection, (2) a dimensionality reduction network that compresses the feature representation into a form suitable for quantum processing, and (3) a hybrid quantum-classical reservoir computing architecture for antigenic distance prediction. Experiments on two datasets covering 1963-2002 (271 strains, 73,441 pairs) and 2003-2025 (888 strains, 788,544 pairs) show that FluQRC outperforms four established baselines across all three evaluation metrics (MAE, RMSE, R2). On the larger and more challenging 2003-2025 dataset, FluQRC achieves MAE = 0.369, RMSE = 0.635, and R2 = 0.900, corresponding to a 20.3% reduction in MAE and a 14.3% reduction in RMSE relative to the strongest baseline, while raising R2 from 0.862 to 0.900. These results demonstrate the scalability and effectiveness of FluQRC for large-scale antigenic distance prediction.

Transcranial direct current stimulation (tDCS) is a promising non-invasive therapy for drug-resistant epilepsy, but the montage optimization pipelines used to personalize it often rely on biophysical head models that capture only the induced electric field and ignore the network dynamics through which stimulation influences epileptogenic activity. We aimed to develop and evaluate a montage optimization framework that incorporates individualized network dynamics. We introduce a framework based on a digital twin, or neurotwin, a personalized replica of a patient's brain that couples a biophysical head model of the stimulation-induced field with a personalized whole-brain model (WBM) of seizure dynamics constrained by structural connectivity and intracranial recordings. We applied it to produce optimal montages in 12 patients with drug-resistant epilepsy, compared them with conventional biophysics-based optimization, and assessed the framework retrospectively against the clinical outcomes of six treated patients. Because the neurotwin optimization suppresses simulated spread by construction, the informative result is the divergence between pipelines: the network-informed and field-based optimizations produced substantially different montages despite delivering comparable inhibitory fields at the epileptogenic zone, and montages with nearly identical field distributions could yield markedly different network outcomes within the model. A minimum-replacement-set analysis traced these differences to field changes at a few propagation-zone and off-target parcels rather than to the overall field. The parcels whose inhibition most reduced spread tended to occupy highly connected (hub) positions in the structural connectome, a robust but partial association indicating that network structure shapes, but does not fully determine, where stimulation is most effective. In the retrospective analysis, model-predicted spread reduction tracked clinical seizure-frequency change in the expected direction, though not significantly (n=6). We hypothesize that accounting for individualized network dynamics may improve seizure control, a hypothesis now being tested prospectively (NCT06334952).